--- title: "3. Integration of rgee with r-spatial ecosystem" author: "Gabriel Carrasco & Antony Barja" date: "2023-01-03 08:28:12" tags: ["rspatial","tmap"] ---
In recent years the R spatial community has had a big impact on big data, especially with the processing of spatial data, both in vector format (points, lines, polygons, etc.) and raster format (satellite images, drones, etc.).
Nowadays, there are many packages within R spatially to work on a key aspect within spatial analysis, from the evaluation of a Moran test to cluster identification, geographic weighted regression, pre and post processing of satellite images or drone images, among others; there is also a great potential within spatial data visualizations, from simple static visualizations, to dynamic and interactive ones, including 3D visualizations.
However, there is a huge gap for processing large volumes of data without high computational costs, and it is precisely in this aspect that the rgee package aims to work and be part of this ecosystem to process large spatial datasets in an accessible and user-friendly way hand in hand with the integration of the spatial universe of R packages.
For this example, we use the ERA5-Land dataset, a reanalysis dataset providing a consistent view of the evolution of land variables over several decades at an enhanced resolution compared to ERA5, more information click here. together MOD13A1 V6 product provides a Vegetation Index, more information click here.
Information: |
- For this demo you need rgee, tmap, raster, viridis,cptcity and mapview packages previously installed. |
library(rgee)
library(tmap)
library(raster)
library(viridis)
library(mapview)
library(cptcity)
ee_Initialize() # Initialize gee inside R with registered account
── rgee 1.1.5 ────────────────────────────────────────────────────── earthengine-api 0.1.334 ──
✔ user: not_defined
✔ Initializing Google Earth Engine: DONE!
✔ Earth Engine account: users/ambarja
───────────────────────────────────────────────────────────────────────────────────────────────
r_temp <- ee$ImageCollection$Dataset$ECMWF_ERA5_LAND_MONTHLY
ee_temp <- r_temp$
select("temperature_2m")$
filterDate("1992-01-01","1992-12-31")$
mean()$
subtract(273.15)
viz_temp <- list(
min = 5,
max = 25,
palette = magma(n = 100)
)
temp <- Map$addLayer(ee_temp, visParams = viz_temp) +
Map$addLegend(visParams = viz_temp,name = "Tem[C°]",position = "bottomleft")
temp
geometry <- ee$Geometry$Rectangle(
coords = c(-180,-90,180,90),
proj = "EPSG:4326",
geodesic = FALSE
)
(minmax <- ee_temp$reduceRegion(
reducer = ee$Reducer$minMax(),
geometry = geometry,
scale = 500*1000)$
getInfo())
$temperature_2m_max
[1] 29.68718
$temperature_2m_min
[1] -53.14849
world_temp <- ee_as_thumbnail(
image = ee_temp, # Image to export
region = geometry, # Region of interest
dimensions = 1024, # output dimension
raster = TRUE, # If FALSE returns a stars object. FALSE by default
vizparams = list(min = -53.14849, max = 29.68718)
)
# Define a projection (EPSG:4326 - WGS 84)
crs(world_temp) <- 4326
world_temp[] <- scales::rescale(
getValues(world_temp),
c(minmax$temperature_2m_min, minmax$temperature_2m_max)
)
# Load world vector (available after load tmap)
data("World")
# Define a Robin type projection
robin <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
# Remove the temperature value of sea
map <- world_temp %>%
mask(World)
tm_shape(map) +
tm_raster(palette = viridis(n = 100),style = "cont")
temp_map <- map %>%
tm_shape(projection = robin,raster.warp = FALSE) +
tm_raster(
title = "Temperature (°C)",
palette = magma(n = 100),
style = "cont"
) +
tm_shape(shp = World) +
tm_borders(col = "black",lwd = 0.7) +
tm_graticules(alpha = 0.8, lwd = 0.5, labels.size = 0.5)
temp_map
tmap_save(temp_map,"temp_map.png",width = 6,height = 3)
# Processing data with rgee
r_ndvi <- ee$ImageCollection$Dataset$MODIS_006_MYD13Q1
ee_ndvi <- r_ndvi$
select("NDVI")$
filterDate("2022-01-01","2022-01-31")$
max()$
multiply(0.0001)
viz_ndvi <- list(
min = -0.5 ,
max = 0.87,
palette = cpt(pal = "grass_ndvi")
)
ndvi <- Map$addLayer(ee_ndvi,visParams = viz_ndvi) +
Map$addLegend(visParams = viz_ndvi,name = "NDVI")
ndvi
temp | ndvi